############################################ ############################################
  ######## REPLICATION FILE FOR: The Value of Dignity Appeals   ########
############################################ ############################################ 


#####################################
####### PRELIMINARIES 
rm(list=ls())
## packages
library(tidyverse)
library(stargazer)
library(texreg)
library(sandwich)
library(lmtest)
library(lmerTest)
library(lme4)
library(MASS)
library(broom.mixed)
library(estimatr)


## set directory to the path that the replication file is in 
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

## create directories in your folder for the figures and tables
dir.create("figures")
dir.create("tables")

## load data
load(file = "repData.RData")
  #loads four datasets: 
    # full_types_i: dataset showing facebook reaction of user  
    # full_main_i: dataset showing facebook any engagement of user
    # comments: dataset with coding of facebook comments
    # pilot: dataset with pilot data on blameworthiness of different groups

####################################
####### FUNCTIONS FOR ANALYSES

graph_tables <- function(dv, data, model){
  # function that takes the data and creates table of regression estimates for the specified model 
  # uses heteroskedasticity robust errors ("HC2")
  # dv: identify the dependent variable of interest
  # data: dataset of interest
  # model: which model you want the regression estiamtes for. Options: 
    # "base": model with treatment variable, controlling for domain (does not consider economic treatment arm)
    # "int": model with treatment interacted with domain (does not consider economic treatment arm)
    # "homeless": model with just treatment (does not consider incarcerated domain)
  
  if(model == "base"){
    
    #if int is the model specified: 
      #relevel the treatment so that info is the reference treatment group
      #filter out those in the economic treatment arm, so that we're comparing info v dignity for incarcerated and homeless domains
      #estimate the effect of treatment controlling for domain
      #get the estimated coefficients and ses 
    
    data$treatment <- relevel(data$treatment, ref = "Info")
    form <- as.formula(paste0(dv, "~ treatment + domain"))
    reg <- lm(form, data = data %>% filter(treatment != "Econ"))
    out <- as_tibble(t(coeftest(reg, vcovHC(reg, type = "HC2"))[c("treatmentDignity"),1:2]))
    out$var <- c("Dignity")
  }else if(model == "int"){
    
    #if int is the model specified: 
      #relevel the treatment so that info is the reference treatment group
      #filter out those in the economic treatment arm, so that we're comparing info v dignity for incarcerated and homeless domains
      #interact domain and treatment to see if there are differences in the effect of dignity across domain
      #get the estimated coefficients and ses 
      
    data$treatment <- relevel(data$treatment, ref = "Info")
    form <- as.formula(paste0(dv, "~ treatment*domain"))
    reg <- lm(form, data = data %>% filter(treatment != "Econ"))
    out <- as_tibble(coeftest(reg, vcovHC(reg, type = "HC2"))[c("treatmentDignity", "treatmentDignity:domainIncarcerated"),1:2])
    out$var <- c("Dignity:Homeless", "Dignity:Incarcerated")
    out$Estimate[out$var == "Dignity:Incarcerated"] <- out$Estimate[out$var == "Dignity:Homeless"] + out$Estimate[out$var == "Dignity:Incarcerated"]
    varOut <- vcovHC(reg, type = "HC2")
    out$`Std. Error`[out$var == "Dignity:Incarcerated"] <-  sqrt(varOut["treatmentDignity", "treatmentDignity"] + 
                                                                   varOut["treatmentDignity:domainIncarcerated", 
                                                                          "treatmentDignity:domainIncarcerated"] + 
                                                                   2*varOut["treatmentDignity", "treatmentDignity:domainIncarcerated"])
  }else if(model == "homeless"){
    
    #if homeless is the model specified: 
      #relevel the treatment so that dignity is the reference treatment group
      #filter out those in the incarcerated domain group, so that just those in the homeless domain group are in the data
      #get the estimated coefficients and ses 
    
    data$treatment <- relevel(data$treatment, ref = "Dignity")
    form <- as.formula(paste0(dv, "~ treatment"))
    reg <- lm(form, data = data %>% filter(domain != "Incarcerated"))
    out <- as_tibble(coeftest(reg, vcovHC(reg, type = "HC2"))[c("treatmentInfo", "treatmentEcon"),1:2])
    out$var <- c("Info", "Econ")
  }
  
  #create confidence intervals, specify the dv, and specify the model
  out$upr <- out$Estimate + 1.96*out$`Std. Error`
  out$lwr <- out$Estimate - 1.96*out$`Std. Error`
  out$fit <- out$Estimate
  out$dv <- dv 
  out$model <- model
  out #output the dataset 
}




##########################################################################
####### MAIN ANALYSES

# this portion of the replication file creates the figures in the main body of the paper

###### FIGURE 2: Effect of dignity on ad engagement

tabs_to_Graph <- lapply(list("pos_reac", "engagement", "anyReac"), function(x){
  #run each model type for every dependent variable of interest: 
    #positive reaction (pos_reac), any engagement (engagement), any reaction (anyReac)
  
  if(x == "engagement"){
    bind_rows(graph_tables(dv = x, data = full_main_i, model = "base"),
              graph_tables(dv = x, data = full_main_i, model = "int"),
              graph_tables(dv = x, data = full_main_i, model = "homeless"))
  }else{
    bind_rows(graph_tables(dv = x, data = full_types_i, model = "base"),
              graph_tables(dv = x, data = full_types_i, model = "int"),
              graph_tables(dv = x, data = full_types_i, model = "homeless"))
  }
  
})

tabs_to_Graph <- do.call("rbind", tabs_to_Graph)

tabs_to_Graph <- tabs_to_Graph %>% mutate(dv_label = case_when(dv == "anyReac" ~ "Reacted to Post",
                                                               dv == "engagement" ~ "Engaged with Post",
                                                               dv == "pos_reac" ~ "Positive Reaction to Post"),
                                          var_label = case_when(var == "Dignity" ~ "Dignity Frame\n(Pooled Results)",
                                                                var == "Econ" ~ "Economic Frame",
                                                                var == "Info" ~ "Information Control",
                                                                var == "Dignity:Homeless" ~ "Dignity Frame\nfor Homeless",
                                                                var == "Dignity:Incarcerated" ~ "Dignity Frame\nfor Incarcerated"))
  #specify the dependent variable and treatment variable labels


#plots the estimated effects and their associated confidence intervals for the 
  #main models - the treatment controlling for the domain and the treatment interacting with domain 
all_plots <- tabs_to_Graph %>% filter(model %in% c("base", "int")) %>% 
  ggplot(aes(x = var_label, y = Estimate, color = dv_label, shape = dv_label)) + 
  geom_point(position =position_dodge(width = 0.5), size = .5) + 
  geom_pointrange(aes(ymin=lwr, ymax=upr), position = position_dodge(width = 0.5)) +
  theme_bw() + 
  geom_hline(yintercept = 0, linetype="dotted", color = "red") +
  labs(title = "Effect of Dignity Frame on Ad Engagement",
       subtitle = "Results from OLS regression",
       y = "Percentage Change in \n User Engagement with Post",
       x = "", color="Dependent Variable", shape = "Dependent Variable") + #coord_flip() + 
  scale_color_manual(values = c("#4097AA", "#2D4057", "#FA9483")) + 
  theme(legend.position = "bottom", 
        title = element_text(size = 28),
        strip.text.x = element_text(size = 20),
        strip.text.y = element_text(size = 20),
        axis.text.y = element_text(size = 24),
        axis.text.x = element_text(size = 24),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 24),
        legend.title = element_text(size = 20),
        legend.text = element_text(size = 20))

ggsave("figures/all_results.png", 
       all_plots, width = 14, height = 6)
  #save the plot


#####################################
####### Appendix
# this portion of the replication file provides code for the 

## ## ## ## ## ## ## ## ## ## ## ## 
## ## ##  Table A1: Summary Statistics 

#creates a summary table for the treatment groups 

sumStats <- full_types_i %>% unite("treat", treatment, domain, sep = ":") %>% #combine frame treatment and domain 
  mutate(`Dignity:Homeless Treatment` = case_when(treat == "Dignity:Homeless" ~ 1, 
                                                  .default = 0),
         `Dignity:Incarcerated Treatment` = case_when(treat == "Dignity:Incarcerated" ~ 1, 
                                                      .default = 0),
         `Econ:Homeless Treatment` = case_when(treat == "Econ:Homeless" ~ 1, 
                                               .default = 0),
         `Info:Homeless Treatment` = case_when(treat == "Info:Homeless" ~ 1, 
                                               .default = 0),
         `Info:Incarcerated Treatment` = case_when(treat == "Info:Incarcerated" ~ 1, 
                                                   .default = 0)) %>% 
  dplyr::select(-campaign_loc, -noReac) %>% #identify variables of interest
  cbind(full_main_i %>% dplyr::select(engagement)) #combine with engagement variable

stargazer(as.data.frame(sumStats),type="latex",title="Summary Statistics", out = "tables/tabA1.tex")
  #creates summary stats for engagement and reactions




## ## ## ## ## ## ## ## ## ## ## ## 
## ## ##  Figure A1: Perceived blameworthiness of different stigmatized groups

# groups of interest
group_names <- c("homeless", "poor", "incarcerated", "disabled", "immigrants", "drugadd")

# Stack the blame variables into long format, count occurrences, and reshape
blame_plot <- pilot %>%
  select_at(all_of(group_names)) %>%  
  drop_na() %>% 
  pivot_longer(cols = everything(),   # Convert to long format
               names_to = "group",
               values_to = "response") %>%
  count(group, response) %>%          # Count occurrences
  pivot_wider(names_from = "group",   # Reshape to wide format
              values_from = "n",
              values_fill = list(n = 0)) #%>%  # Fill missing counts with 0
 # rename_with(~ group_names, -response) # Rename columns

blame_plotl <- blame_plot %>% 
  pivot_longer(-response) %>% #make the data in long format
  filter(! response == "Missing") %>% #get rid of missing data
  mutate(Labels2 = case_when( #create ordered labels for plotting
    response == "Neutral" ~ 3,
    response == "Not responsible at all" ~ 1,
    response == "Not very responsible" ~ 2,
    response == "Somewhat responsible" ~ 4,
    response == "Very responsible" ~ 5
  ))

#create bar graph of blameworthiness plot
(ggplot(blame_plotl, aes(x = Labels2, y = value, fill = name)) + geom_bar(stat = "identity", position = "dodge") + 
  facet_wrap(~ name) + theme_bw() + 
  labs(fill = "Group", y = "Frequency", x = "Blameworthiness (1 = Least Blame)") + 
  theme(legend.position = "top") ) %>% 
ggsave("figures/blame.png",., width = 10, height = 6)



## ## ## ## ## ## ## ## ## ## ## ## 
## ## ##  Table A2: Table of results
#creates the table version of the results from figure 2 in main analysis

#relevel the treatment variable so information is the reference group
full_types_i$treatment <- relevel(factor(full_types_i$treatment, ordered = FALSE), ref = "Info")
full_main_i$treatment <- relevel(factor(full_main_i$treatment, ordered = FALSE), ref = "Info")

## base model - dignity frame v info, controlling for domain (excludes econ frame)
# (POSITIVE REACTION)
posBase <- lm_robust(pos_reac ~ treatment + domain , full_types_i[full_types_i$treatment != "Econ",], se_type = "HC2")
# (ANY REACTION)
reacBase <-lm_robust(anyReac ~ treatment +  domain , full_types_i[full_types_i$treatment != "Econ",], se_type = "HC2")
# (ANY ENGAGEMENT)
engageBase <- lm_robust(engagement ~ treatment + domain , full_main_i[full_main_i$treatment != "Econ",], se_type = "HC2")

## int model - dignity frame v info, interacted with domain (excludes econ frame)
# (POSITIVE REACTION)
posBaseI <- lm_robust(pos_reac ~ treatment*domain , full_types_i[full_types_i$treatment != "Econ",])
# (ANY REACTION)
reacBaseI <-lm_robust(anyReac ~ treatment*domain , full_types_i[full_types_i$treatment != "Econ",])
# (ANY ENGAGEMENT)
engBaseI <- lm_robust(engagement ~ treatment*domain , full_main_i[full_main_i$treatment != "Econ",])


### Homeless -> Dignity v. Control v. Econ
#relevels treatment so dignitiy is the reference category
full_types_i$treatment <- relevel(factor(full_types_i$treatment, ordered = FALSE), ref = "Dignity")
full_main_i$treatment <- relevel(factor(full_main_i$treatment, ordered = FALSE), ref = "Dignity")

# (POSITIVE REACTION)
posBaseH <- lm_robust(pos_reac ~ treatment , full_types_i[full_types_i$domain != "Incarcerated",])
# (ANY REACTION)
reacBaseH <-lm_robust(anyReac ~ treatment , full_types_i[full_types_i$domain != "Incarcerated",])
# (ANY ENGAGEMENT)
engBaseH <- lm_robust(engagement ~ treatment , full_main_i[full_main_i$domain != "Incarcerated",])


### combine results into one table
tabA2 <- texreg(l = list(posBase, posBaseI, posBaseH, reacBase, reacBaseI, reacBaseH, engageBase, engBaseI, engBaseH), 
       include.ci = FALSE, digits = 4,
       stars = c(0.001, 0.01, 0.05, 0.1),
       custom.coef.names = c("Intercept", "Dignity Frame (Base = Information)", 
                             "Vulnerable Group Incarcerated",
                             "Dignity Frame:Incarcerated",
                             "Information Frame (Base = Dignity)",
                             "Economic Frame (Base = Dignity)"),
       custom.model.names = c("(1)", "(2)","(3)", "(4)","(5)", "(6)", "(7)", "(8)", "(9)"),
       custom.gof.rows = list("Domain" = c("Both", "Both", "Homeless Only" ,
                                           "Both", "Both", "Homeless Only",
                                           "Both", "Both", "Homeless Only")),
       custom.header = list("Positive Reaction" = 1:3, "Any Reaction" = 4:6, "Any Engagement" = 7:9))

writeLines(tabA2, "tables/tabA2.tex")


## ## ## ## ## ## ## ## ## ## ## ## 
## ## ## Section A7.1 

#calculates adjusted p-values for BH corrections for the three sets of models. (Referenced in the paper)
all_pesTMP <- lapply(as.list(c("posBase", "reacBase", "engageBase",
                               "posBaseI", "reacBaseI", "engBaseI",
                               "posBaseH", "reacBaseH", "engBaseH")), 
                     function(x) 
                       
                     {
                       out <- broom::tidy(get(x)) %>% 
                         filter(! term == "(Intercept)") %>% 
                         mutate(adjust = p.adjust(p.value, method = "fdr", n = 3)) %>% 
                         mutate(reference_name = x)
                     }
)

all_pes <- do.call(rbind, all_pesTMP)


## ## ## ## ## ## ## ## ## ## ## ## 
## ## ##  Table A3: Multi-level models

#estimates multi-level models accounting for the structure of the treatment assignment 

full_main_i <- within(full_main_i, treatment <- relevel(treatment, ref = "Info"))
full_types_i <- within(full_types_i, treatment <- relevel(treatment, ref = "Info"))

pooled_noecon_mlm_eng <- lmer(engagement ~ treatment + domain + (1 | campaign_loc), 
                              data = full_main_i[full_main_i$treatment != "Econ",])
class(pooled_noecon_mlm_eng) <- "lmerMod"
pooled_noecon_mlm_anyr <- lmer(anyReac ~ treatment + domain + (1 | campaign_loc), 
                               data = full_types_i[full_types_i$treatment != "Econ",])
class(pooled_noecon_mlm_anyr) <- "lmerMod"
pooled_noecon_mlm_pos <- lmer(pos_reac ~ treatment + domain + (1 | campaign_loc), 
                              data = full_types_i[full_types_i$treatment != "Econ",])
class(pooled_noecon_mlm_pos) <- "lmerMod"

#### pooled and interactive effects without economic 
pooled_noecon_mlm_int_eng <- lmer(engagement ~ treatment*domain + (1 | campaign_loc), 
                                  data = full_main_i[full_main_i$treatment != "Econ",])
class(pooled_noecon_mlm_int_eng) <- "lmerMod"
pooled_noecon_mlm_int_anyr <-lmer(anyReac ~ treatment*domain + (1 | campaign_loc), 
                                  data = full_types_i[full_types_i$treatment != "Econ",])
class(pooled_noecon_mlm_int_anyr) <- "lmerMod"
pooled_noecon_mlm_int_pos <- lmer(pos_reac ~ treatment*domain + (1 | campaign_loc), 
                                  data = full_types_i[full_types_i$treatment != "Econ",])
class(pooled_noecon_mlm_int_pos) <- "lmerMod"

# focusing only on homelessness
full_main_i <- within(full_main_i, treatment <- relevel(treatment, ref = "Dignity"))
full_types_i <- within(full_types_i, treatment <- relevel(treatment, ref = "Dignity"))

dignity_mlm_eng <- lmer(engagement ~ treatment + (1 | campaign_loc), 
                        data = full_main_i[full_main_i$domain != "Incarcerated",])
class(dignity_mlm_eng) <- "lmerMod"


dignity_mlm_anyr <- lmer(anyReac ~ treatment + (1 | campaign_loc), 
                         data = full_types_i[full_main_i$domain != "Incarcerated",])
class(dignity_mlm_anyr) <- "lmerMod"


dignity_mlm_pos <- lmer(pos_reac ~ treatment + (1 | campaign_loc), 
                        data = full_types_i[full_main_i$domain != "Incarcerated",])
class(dignity_mlm_pos) <- "lmerMod"

mlmMods <- as.list(c(dignity_mlm_pos, pooled_noecon_mlm_pos, pooled_noecon_mlm_int_pos, 
                  dignity_mlm_anyr, pooled_noecon_mlm_anyr, pooled_noecon_mlm_int_anyr,
                  dignity_mlm_eng, pooled_noecon_mlm_eng, pooled_noecon_mlm_int_eng))
  #create a list of the mlm model objects

stargazer(mlmMods,
          dep.var.labels = c("Positive Reaction", "Any Reaction", "Any Engagement"), 
          out = "tables/tabA3.tex")
  #create a table of the mlm models



## ## ## ## ## ## ## ## ## ## ## ## 
## ## ##  Figure A2: Effect of dignity frame versus info and econ

#plots the effect for just the homeless domain looking at the difference between dignity v info and dignity v econ framing

baseECON <- tabs_to_Graph %>% filter(model == "homeless") %>% #runs regression and specifies that the homeless model is of interest
  mutate(Estimate = -1*Estimate, #flips the direction so that it is the effect of dignity v info or econ treatment
         lwr = -1*lwr,
         upr = -1*upr,
         var_label = case_when(var == "Econ" ~ "Dignity v. Economic Frame",
                               var == "Info" ~ "Dignity v. Information Control")) %>% 
  ggplot(aes(x = var_label, y = Estimate, color = dv_label)) + 
  geom_point(position =position_dodge(width = 0.5), size = .5) + 
  geom_pointrange(aes(ymin=lwr, ymax=upr), position = position_dodge(width = 0.5)) +
  theme_bw() + 
  geom_hline(yintercept = 0, linetype="dotted", color = "red") +
  labs(title = "Effect of Dignity Frame versus Control and Economic Frame",
       subtitle = "OLS regression",
       y = "Percentage Change in \n User Engagement with Post",
       x = "", color="Dependent Variable") + #coord_flip() + 
  scale_color_manual(values = c("#4097AA", "#2D4057", "#FA9483")) + 
  theme(legend.position = "top", 
        title = element_text(size = 20),
        strip.text.x = element_text(size = 15),
        strip.text.y = element_text(size = 15),
        axis.text.y = element_text(size = 14),
        axis.text.x = element_text(size = 14),
        axis.title.x = element_text(size = 17),
        axis.title.y = element_text(size = 17),
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 14))

ggsave("figures/baseHOMELESS.png", baseECON, width = 12, height = 6)


## ## ## ## ## ## ## ## ## ## ## ## 
## ## ##  Table A4: summary stats for coding of content in facebook comments
stargazer(as.data.frame(comments %>% dplyr::select(Relevant:Blameworthy)), out = "tables/tabA4.tex")
  #summary statistic table for the comments dataset


## ## ## ## ## ## ## ## ## ## ## ## 
## ## ##  Table A5: comment effects
#estimates the same models for comments as the outcome

# base relationships for comments (effect of frame (dignity v info) controlling for domain)
comments$Treatment <- relevel(factor(comments$Treatment, ordered = FALSE), ref = "Information")
base_pr <- lm(Positive_Respect ~ Treatment + Topic, data = comments[comments$Treatment != "Econ",])
base_nd <- lm(Negative_Disrespect ~ Treatment + Topic, data = comments[comments$Treatment != "Econ",])
base_b <- lm(Blameworthy ~ Treatment + Topic, data = comments[comments$Treatment != "Econ",])


# interaction term for comments (effect of frame (dignity v info) interacting domain)
int_pr <- lm(Positive_Respect ~ Treatment*Topic, data = comments[comments$Treatment != "Econ",])
int_nd <- lm(Negative_Disrespect ~ Treatment*Topic, data = comments[comments$Treatment != "Econ",])
int_b <- lm(Blameworthy ~ Treatment*Topic, data = comments[comments$Treatment != "Econ",])


# homelessness (effect of dignity (v info and v econ) for just homelessness)
comments$Treatment <- relevel(comments$Treatment, ref = "Dignity")

h_pr <- lm(Positive_Respect ~ Treatment, data = comments[comments$Topic != "Incarcerated",])
h_nd <- lm(Negative_Disrespect ~ Treatment, data = comments[comments$Topic != "Incarcerated",])
h_b <- lm(Blameworthy ~ Treatment, data = comments[comments$Topic != "Incarcerated",])

# create a table with these results
tabA5 <- texreg(l = list(base_pr, int_pr, int_pr, base_nd, int_nd, h_nd, base_b, int_b, h_b), 
       include.ci = FALSE, digits = 4,
       stars = c(0.001, 0.01, 0.05, 0.1),
       custom.coef.names = c("Intercept", "Dignity Frame (Base = Information)", 
                             "Vulnerable Group Incarcerated",
                             "Dignity Frame:Incarcerated",
                             "Information Frame (Base = Dignity)",
                             "Economic Frame (Base = Dignity)"),
       custom.model.names = c("(1)", "(2)","(3)", "(4)","(5)", "(6)", "(7)", "(8)", "(9)"),
       custom.gof.rows = list("Domain" = c("Both", "Both", "Homeless Only" ,
                                           "Both", "Both", "Homeless Only",
                                           "Both", "Both", "Homeless Only")),
       custom.header = list("Positive or Respect" = 1:3, "Negative or Disrespect" = 4:6, "Blameworthy" = 7:9))

writeLines(tabA5, "tables/tabA5.tex")
